
rm(list=ls(all=TRUE))

setwd("replication_archive/MonteCarlo/output")

library(MASS)

set.seed(123)

beta0 <- 1
beta1 <- 2
beta2 <- -3
beta3 <- 1

data <- NULL

for(i in 1:1000){

	bla <- mvrnorm(1000, mu=c(3, 1), Sigma=matrix(c(2, 0.3, 0.3, 1), nrow=2))
	t1 <- bla[,1]
	w1 <- bla[,2]
	c1 <- rbinom(1000, 1, 0.2)
	e1 <- rnorm(1000, 0, 1)

	y <- beta0 + beta1*t1 + beta2*c1 + beta3*w1 + e1

	m_t <- lm(y ~ t1 + c1 + w1)

	x1 <- ((1-c1) * t1) + (c1*(t1*2))

	m_x <- lm(y ~ x1 + c1 + w1)

	data <- rbind(data, coef(m_x))
}

data <- as.data.frame(data)

colnames(data) <- c("delta0", "delta1", "delta2", "delta3")


quartz(type="pdf", width=4, height=4, file="mc_delta1.pdf")
par(mar = c(2,2,0.1,0.1), mgp=c(2,0.5,0), family="CMU Serif")
plot(density(data$delta1), lwd=5, xlab="", main="", ylab="", xlim=c(min(density(data$delta1)$x), beta1))
abline(v=beta1, lwd=5, col="red")
dev.off()

quartz(type="pdf", width=4, height=4, file="mc_delta2.pdf")
par(mar = c(2,2,0.1,0.1), mgp=c(2,0.5,0), family="CMU Serif")
plot(density(data$delta2), lwd=5, xlab="", main="", ylab="", xlim=c(min(density(data$delta2)$x), beta2))
abline(v=beta2, lwd=5, col="red")
dev.off()

quartz(type="pdf", width=4, height=4, file="mc_delta3.pdf")
par(mar = c(2,2,0.1,0.1), mgp=c(2,0.5,0), family="CMU Serif")
plot(density(data$delta3), lwd=5, xlab="", main="", ylab="", xlim=c(min(density(data$delta3)$x), max(density(data$delta3)$x)))
abline(v=beta3, lwd=5, col="red")
dev.off()

